Introduction

This exercise will build on the skills acquired in the Mapping Basics and Species Distribution Modeling exercises to perform Home Range Analysis and map the results. We will start by examining the structure of data used in home range analyses, we will look at ways to perform QA/QC checks on the location data, and finally various ways to visualize the data. One of the new skills we will discuss in this exercise is to create functions and loops to replicate analyses.

Packages used in this exercise

There are several specialty packages that will be used in this exercise due to the specific nature of the analyses. Some of these packages you will need to install while several others we have used in previous exercises and should already be installed.

To begin, we will install the following:

library("adehabitatHR")
library("data.table")
library("ggfortify")
library("grid")
library("move")
library("moveVis")
library("pbapply")
library("plotly")
library("rgdal")
library("sp")
library("tidyverse")
library("viridis")
library('ggmap')
library("OpenStreetMap")
library("gpclib")
library("maptools")

Dataset

Importing dataset with read.csv:

data <- read.csv("./Data/lowaSCH.csv")

This data comes from my thesis research from this past spring/summer. This dataset contains only my study site on the street behind Dr. Schillers’ property, which contains GPS points for 3 individual male Louisiana Waterthrushes.

qaqc_plot <- ggplot() + geom_point(data=data, 
                                   aes(location.long,location.lat,
                                       color=individual.local.identifier)) +
                        labs(x="Longitude", y="Latitude") +
                        guides(color=guide_legend("Individual"))

ggplotly(qaqc_plot)

With plotly, similar to leaflet, we have the ability to examine the spread of the data points and additional information from various columns. From this plot we can see that there are three individuals; m1, m2, and m3. Already you can see how their territories follow pretty tightly along the stream and tributary at this site.

While we could continue with the current dataset, any analysis would calculate home range for the entire population rather than the individual.

Creating a function and using the lapply command to apply a function over a list or vector dataset. Specifically, this function will take the original dataset, split it into separate files based on the individual identifier, and create new *.csv files using the identifier as the filename.

lapply(split(data, data$individual.local.identifier), 
       function(x)write.csv(x, file = paste(x$individual.local.identifier[1],".csv", sep = ""), row.names = FALSE))

files <- list.files(path = ".", pattern = "[m]+[0-9]+", full.names = TRUE)

In the list.files command above, path = "." informs the locations, in this case the root directory, pattern = describes the way the files are named, in this case “m” (for male) followed by a number between 0-9, and full.names describes how the files will be listed.

Analysis

Imagery

Although my data contained only longitude/latitude values, I converted them to UTM and back again to longitude/latitude coordinates to practice using the code if I ever receive data that only has UTM values.

utm_points <- cbind(data$utm.easting, data$utm.northing)
utm_locations <- SpatialPoints(utm_points, 
                 proj4string=CRS("+proj=utm +zone=16 +datum=WGS84"))
proj_lat.lon <- as.data.frame(spTransform(
                utm_locations, CRS("+proj=longlat +datum=WGS84")))
colnames(proj_lat.lon) <- c("x","y")
raster <- openmap(c(max(proj_lat.lon$y)+0.001, min(proj_lat.lon$x)-0.001), 
                  c(min(proj_lat.lon$y)-0.001, max(proj_lat.lon$x)+0.001), 
                  type = "bing")
raster_utm <- openproj(raster, 
              projection = "+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs")

In the script above, utm_point is an x,y derived from the primary dataset, utm_locations set the projection to UTM Zone 16, proj_lat.lon converted the UTM points to longitude/latitude, raster uses the min/max x,y data to create a bounding box to retrieve the aerial imagery, and raster_utm reprojected the imagery back to UTM Zone 16 for Tennessee. Next, autoplot.OpenStreetMap will display the raster image file with the UTM locations as an overlay.

autoplot.OpenStreetMap(raster_utm, expand = TRUE) + theme_bw() +
  theme(legend.position="bottom") +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_point(data=data, aes(utm.easting,utm.northing,
             color=individual.local.identifier), size = 3, alpha = 0.8) +
  theme(axis.title = element_text(face="bold")) + labs(x="Easting",
        y="Northing") + guides(color=guide_legend("Identifier"))

Home Range Analysis

In the section above we use the lapply command to loop a function used to separate the original dataset into individual files. This is a useful tool, however, when the function loops through dozens or even hundreds of files, the process can take a long period of time to complete. Using the pblapply command adds a progress bar (i.e. pb) to the process which provides an estimated time for completion of the function. We will use a similar process to the one above, using the pblapply command to run MCP analysis on the individual *.csv files.

Minimum Convex Polygon

A description of the steps within the following code will be discussed following the output. The process works by establishing the function and then running pblapply referring back to the files we created in the dataset portion of this exercise.

I had some issues getting “pblapply” to run, so I had to install the package “gpclib” and specify type=“source” in order to bypass whatever error I was getting.

I also changed the point and polygon colors in order to see them better on the dark background (especially since I couldn’t get it to zoom in any closer to the points).

if (!require(gpclib)) install.packages("gpclib", type="source")
gpclibPermit()

## [1] TRUE

mcp_raster <- function(filename){
  data <- read.csv(file = filename)
  x <- as.data.frame(data$utm.easting)
  y <- as.data.frame(data$utm.northing)
  xy <- c(x,y)
  data.proj <- SpatialPointsDataFrame(xy,data, proj4string = CRS("+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs"))
  xy <- SpatialPoints(data.proj@coords)
  mcp.out <- mcp(xy, percent=100, unout="ha")
  mcp.points <- cbind((data.frame(xy)),data$individual.local.identifier)
  colnames(mcp.points) <- c("x","y", "identifier")
  mcp.poly <- fortify(mcp.out, region = "id")
  units <- grid.text(paste(round(mcp.out@data$area,2),"ha"), x=0.85,  y=0.95,
                     gp=gpar(fontface=4, col="white", cex=0.9), draw = FALSE)
  mcp.plot <- autoplot.OpenStreetMap(raster_utm, expand = TRUE) + theme_bw() + theme(legend.position="none") +
    theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
    geom_polygon(data=mcp.poly, aes(x=mcp.poly$long, y=mcp.poly$lat), color="#ff5e00",fill="#c34800", alpha=0.7) +
    geom_point(data=mcp.points, aes(x=x, y=y), color = "#ff9b00", size=3) + 
    labs(x="Easting (m)", y="Northing (m)", title=mcp.points$identifier) +
    theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5)) + 
    annotation_custom(units)
  mcp.plot
}

pblapply(files, mcp_raster)

## [[1]]

## 
## [[2]]

## 
## [[3]]

Kernel-Density Estimation

The same design above can be used to create KDE plots for each individual in the dataset. Because the process is essentially the same, only portions of the function that were not previously described will be discussed.

kde_raster <- function(filename){
  data <- read.csv(file = filename)
  x <- as.data.frame(data$utm.easting)
  y <- as.data.frame(data$utm.northing)
  xy <- c(x,y)
  data.proj <- SpatialPointsDataFrame(xy,data, proj4string = CRS("+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs"))
  xy <- SpatialPoints(data.proj@coords)
  kde<-kernelUD(xy, h="href", kern="bivnorm", grid=100)
  ver <- getverticeshr(kde, 30)
  kde.points <- cbind((data.frame(data.proj@coords)),data$individual.local.identifier)
  colnames(kde.points) <- c("x","y","identifier")
  kde.poly <- fortify(ver, region = "id")
  units <- grid.text(paste(round(ver$area,2)," ha"), x=0.85,  y=0.95,
                     gp=gpar(fontface=4, col="white", cex=0.9), draw = FALSE)
  kde.plot <- autoplot.OpenStreetMap(raster_utm, expand = TRUE) + theme_bw() + theme(legend.position="none") +
    theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
    geom_polygon(data=kde.poly, aes(x=kde.poly$long, y=kde.poly$lat), color="#0d9be8", fill="#0d9be8", alpha = 0.8) +
    geom_point(data=kde.points, aes(x=x, y=y), color="#6ee1f0", size=4) +
    labs(x="Easting (m)", y="Northing (m)", title=kde.points$identifier) +
    theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5)) + 
    annotation_custom(units)
  kde.plot
}

pblapply(files, kde_raster)

## [[1]]

## 
## [[2]]

## 
## [[3]]

I kept getting errors with this map, saying that “The grid is too small to allow the estimation of home-range. You should rerun kernelUD with a larger extent parameter”.

After trying all kinds of code fixes, what finally got this to run was changing the “getverticeshr(kde, 95)” from 95 to 30. I’m not sure why this is the magic number that worked, but it did give me some home range estimates as you can see in the maps above.

Brownian Bridge Movement

Although this appears to be working backwards, in the previous two examples you have seen how to create a function that can be looped to run the analysis for a list fo files. For this analysis, we will create a conventional script for a single individual. Portions of the script that were not previously described will be discussed following the analysis.

#LOWAsch <- read.csv("./Data/lowaSCH.csv")

#date <- as.POSIXct(strptime(as.character(LOWAsch$timestamp),"%Y-%m-%d %H:%M:%S", tz="US"))
#LOWAsch$date <- date
#LOWAsch.reloc <- cbind.data.frame(LOWAsch$utm.easting, LOWAsch$utm.northing,
#                                as.vector(LOWAsch$individual.local.identifier),
#                                as.POSIXct(date))
#colnames(LOWAsch.reloc) <- c("x","y","id","date")
#trajectory <- as.ltraj(LOWAsch.reloc, date=date, id="LOWAsch")
#sig1 <- liker(trajectory, sig2 = 58, rangesig1 = c(0, 5), plotit = FALSE)
#lowa.traj <- kernelbb(trajectory, sig1 = .7908, sig2 = 58, grid = 100)
#bb_ver <- getverticeshr(lowa.traj, 95)
#bb_poly <- fortify(bb_ver, region = "id", 
#                   proj4string = CRS("+proj=utm +zone=16+
#                                     datum=WGS84 +units=m +no_defs"))
#colnames(bb_poly) <- c("x","y","order","hole","piece","id","group")
#bb_image <- crop(lowa.traj, bb_ver, 
#                 proj4string = CRS("+proj=utm +zone=16 +
#                                   datum=WGS84 +units=m +no_defs"))
#bb_units <- grid.text(paste(round(bb_ver$area,2)," ha"), x=0.85,  y=0.95,
#                   gp=gpar(fontface=4, col="white", cex=0.9), draw = FALSE)
#bb.plot <- autoplot.OpenStreetMap(raster_utm, expand = TRUE) + theme_bw() + theme(legend.position="none") +
#  theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
#  geom_tile(data=bb_image, 
#            aes(x=bb_image@coords[,1], y=bb_image@coords[,2],
#            fill = bb_image@data$ud)) +
#  geom_polygon(data=bb_poly, aes(x=x, y=y, group = group), color = "black", fill = NA) +
#  scale_fill_viridis_c(option = "inferno") + annotation_custom(bb_units) +
#  labs(x="Easting (m)", y="Northing (m)", title="LOWAsch") +
#  theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5))
#bb.plot

Issue w/ dates

Animate Trajectory Data

The final analysis we will perform in this exercise is a visual animation of the relocation data. Using the data from above you can plot(trajectory) and see a plot of the relocation information for an individual. However, we can use the move and moveVis packages to create an animation of the relocations.

We need to begin by creating a move object containing relocations (x,y), time and date information (time), a projection string, individual identifier (animal), and sensor type (sensor). See the description for move() in the help menu for optional information.

#lowa_move <- move(x=LOWAsch$location.long, 
#lowa2_move<-move(df$location.long, df$location.lat,
#             y=LOWAsch$location.lat, 
 #            time=as.POSIXct(LOWAsch$timestamp, 
  #                           format="%m/%d/%Y %H:%M",tz="UTC"), 
             #time=as.POSIXct(df$timestamp,format="%m/%d/%Y %H:%M",tz="UTC"),
   #            proj=CRS("+init=epsg:32615"))
    #         proj=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
     #        data=LOWAsch
      #       animal=LOWAsch$individual.local.identifier 
       #      sensor=LOWAsch$sensor.type

#movement <- align_move(lowa_move, res = "max", digit = 0, unit = "secs")

With the data now on a uniform time scale we can create the frames and animation for the relocations. For this step I will use a basemap from MapBox which require token access through the use of their API. To do this you need to register with MapBox and create an access token. Then create a .Renviron file in your project folder. Copy the token information from MapBox and create an object in the .Renviron file such as map_token = 'paste token here' and add map_token = Sys.getenv('map_token') to the script below. However using the get_maptypes() script you can see there are various map services and map types that can be used. A simple output would be to use map_service = 'osm' (OpenStreetMaps) and map_type = 'topographic' or other map types available by viewing get_maptypes('owm'). when using a basemap without token access the map_token option can be removed from the script below.

#frames <- frames_spatial(movement, path_colors = "red",
 #                        map_service = "osm", 
  #                       map_type = "topographic",
   #                      alpha = 0.5) %>% 
#  add_labels(x = "Longitude", y = "Latitude") %>%
 # add_northarrow() %>% 
  #add_timestamps(movement, type = "label") %>%
#  add_scalebar(distance = 2) %>% 
 # add_progress()

#animate_frames(frames, fps = 5, overwrite = TRUE,
 #              out_file = "./moveVis-2021lowa.gif")

Home Range Analysis
Advanced Data Analytics

Jenna Atma

Introduction

Packages used in this exercise

Dataset

Analysis

Imagery

Home Range Analysis

Minimum Convex Polygon

Kernel-Density Estimation

Brownian Bridge Movement

Animate Trajectory Data

Home Range Analysis Advanced Data Analytics

Jenna Atma

Introduction

Packages used in this exercise

Dataset

Analysis

Imagery

Home Range Analysis

Minimum Convex Polygon

Kernel-Density Estimation

Brownian Bridge Movement

Animate Trajectory Data

Home Range Analysis
Advanced Data Analytics